﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "pqg2")]
    public static unsafe double pqg2(double a, double b, double eps, IntPtr s_x_ya_ptr, IntPtr f_xy_ptr)
    {
        s_x_ya = Marshal.GetDelegateForFunctionPointer<delegatefunc_x_ya>(s_x_ya_ptr);
        f_xy = Marshal.GetDelegateForFunctionPointer<delegatefunc_xy>(f_xy_ptr);

        return pqg2(a, b, eps);
    }

    /// <summary>
    /// 计算二重积分的连分式
    /// 固定一个x，用连分式法计算一个对y的积分近似值。
    /// </summary>
    /// <param name="x"></param>
    /// <param name="eps"></param>
    /// <returns></returns>
    private static unsafe double pqg1(double x, double eps = 1.0E-5)
    {
        int m, n, k, j, flag;
        double h0, g0, s0, s1 = 0.0, d, yy;
        double* h = stackalloc double[10];
        double* g = stackalloc double[10];
        double* b = stackalloc double[10];
        double* y = stackalloc double[2];

        m = 0;
        n = 1;
        // 计算上下限y[1]与y[0]
        s_x_ya(x, y);
        h0 = y[1] - y[0];
        flag = 0;
        // 梯形公式计算初值
        g0 = h0 * (f_xy(x, y[0]) + f_xy(x, y[1])) / 2.0;
        while ((m < 10) && (flag == 0))
        {
            m = m + 1;
            h[0] = h0;
            g[0] = g0;
            // 计算b[0]
            b[0] = g[0];
            j = 1;
            s1 = g[0];
            while (j <= 7)
            {
                d = 0.0;
                for (k = 0; k <= n - 1; k++)
                {
                    yy = y[0] + (k + 0.5) * h[j - 1];
                    d = d + f_xy(x, yy);
                }
                // 变步长梯形求积法计算新近似值g[j]
                g[j] = (g[j - 1] + h[j - 1] * d) / 2.0;
                h[j] = h[j - 1] / 2.0;
                n = 2 * n;
                // 计算b[j]
                funpqj(h, g, b, j);
                s0 = s1;
                // 连分式法计算积分近似值s1
                s1 = funpqv(h, b, j, 0.0);
                if (Math.Abs(s1 - s0) >= eps) j = j + 1;
                else j = 10;
            }
            h0 = h[j - 1]; g0 = g[j - 1];
            if (j == 10) flag = 1;
        }
        return (s1);
    }

    /// <summary>
    /// 计算二重积分的连分式法
    /// s(double,double[])计算上下限的函数名。
    /// f(double,double)计算被积函数值f(x,y)的函数名。
    /// </summary>
    /// <param name="a">a积分下限。</param>
    /// <param name="b">b积分精度要求。积分上限。要求b>a。</param>
    /// <param name="eps"></param>
    /// <returns>函数返回积分值。</returns>
    public static unsafe double pqg2(double a, double b, double eps = 1.0E-5)
    {
        int k, j, m, flag, n;
        double h0, g0, d, s0, s1, x;
        double* h = stackalloc double[10];
        double* g = stackalloc double[10];
        double* bb = stackalloc double[10];

        m = 0; n = 1;
        h0 = b - a; flag = 0;
        // 固定 x = a
        s0 = pqg1(a, eps);
        // 固定 x = b
        s1 = pqg1(b, eps);
        // 梯形公式计算初值
        g0 = h0 * (s1 + s0) / 2.0;
        while ((m < 10) && (flag == 0))
        {
            m = m + 1;
            h[0] = h0;
            g[0] = g0;
            // 计算b[0]
            bb[0] = g[0];
            j = 1;
            s1 = g[0];
            while (j <= 7)
            {
                d = 0.0;
                for (k = 0; k <= n - 1; k++)
                {
                    // 固定一个x
                    x = a + (k + 0.5) * h[j - 1];
                    d = d + pqg1(x, eps);
                }
                // 变步长梯形求积法计算新近似值g[j]
                g[j] = (g[j - 1] + h[j - 1] * d) / 2.0;
                h[j] = h[j - 1] / 2.0;
                n = 2 * n;
                // 计算b[j]
                funpqj(h, g, bb, j);
                s0 = s1;
                // 连分式法计算积分近似值s1
                s1 = funpqv(h, bb, j, 0.0);
                if (Math.Abs(s1 - s0) >= eps) j = j + 1;
                else j = 10;
            }
            h0 = h[j - 1];
            g0 = g[j - 1];
            if (j == 10)
            {
                flag = 1;
            }
        }
        return (s1);
    }

    /*
    // 计算二重积分的连分式法例
      int main()
      { 
          double a,b,eps,s,pqg2f(double,double);
          void  pqg2s(double,double []);
          a=0.0; b=1.0; eps=0.00001;
          s=pqg2(a,b,eps,pqg2s,pqg2f);
          cout <<"s = " <<s <<endl;
          return 0;
      }
    // 计算上下限y1(x)与y0(x)
      void pqg2s(double x, double y[2])
      { 
          y[1]=sqrt(1.0-x*x);
          y[0]=-y[1];
          return;
      }
    // 计算被积函数值f(x,y)
      double pqg2f(double x, double y)
      { 
          double z;
          z=exp(x*x+y*y);
          return(z);
      }
    */
}

